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ABSTRACT 

A  vector  space  approach  to  image  reconstruction  is  derived 
and  introduced.  The  continuous-domain  image  is  assumed 
to  belong  to  a  reproducing  kernel  Hilbert  space  and  the 
sampling  process  is  shown  to  correspond  to  an  appropri¬ 
ate  orthogonal  projection.  The  values  at  the  interpolating 
grid  are  shown  to  correspond  to  a  set  of  inner  product 
calculations,  giving  rise  to  a  minimax  solution  for  an  £2 
approximation  problem.  A  tight  upper  bound  on  the  ensued 
error  is  then  derived  and  demonstrated.  Examples  of  image 
resizing  show  that  the  proposed  method  yields  better  results 
than  presently  available  methods,  including  the  cubic  B- 
spline  method,  in  terms  of  SNR. 


mation  error.  Also,  the  piecewise-polynomial  model  is  gen¬ 
eralized  by  considering  2D  Sobolev  signals  of  finite  support. 
Such  an  assumption  on  the  continuous-domain  image  has 
been  made  in  several  image  processing  algorithms  [7], 
[8]  and  it  is  adopted  here  as  well.  Within  this  setting, 
the  ideal  sampling  process  is  shown  to  correspond  to  a 
set  of  inner  product  calculations;  the  same  holds  for  the 
interpolated  values.  This  signal  representation  interpretation 
to  the  sampling  and  to  the  interpolation  processes  gives  rise 
to  an  alternative  interpolation  approach.  A  detailed  analysis 
is  given  for  the  ID  case,  followed  by  2D  scaling  examples. 

II.  SAMPLING  OF  SMOOTH  SIGNALS 


I.  INTRODUCTION 


Image  interpolation  is  a  fundamental  task  in  image  process¬ 
ing  applications.  Such  applications  include  rotation,  trans¬ 
lation,  resizing  and  derivative  evaluation  to  name  a  few. 
The  underlying  idea  in  current  interpolation  methods  cor¬ 
responds  to  regularity  properties  that  are  assumed  on 
the  continuous-domain  image.  For  example,  a  piecewise- 
polynomial  model  is  often  used,  implying  that  the  original 
continuous-domain  image  is  smooth,  up  to  a  certain  degree. 
Several  interpolation  methods  such  as  nearest-neighbor, 
linear,  Dodgson,  Keys,  Schaum,  B-spline  of  higher  orders, 
Meijering  and  o-MOMS  assume  such  a  model  [1]— [5]. 
Furthermore,  for  a  sufficiently  smooth  input  signal,  these 
very  methods  comply  with  the  following  upper  bound  on 
the  approximation  error 


x  —  x 


l2  - 


<C  -TL  • 


rW 


as  T  — >  0. 


(1) 


Here,  x  is  the  original  continuous-domain  signal,  x  is  the 
interpolated  signal  and  T  is  the  sampling  interval.  In  such 
a  formulation,  the  parameters  L  and  C  are  the  approxima¬ 
tion  order  and  the  proportional  constant,  respectively;  they 
provide  a  means  for  comparing  the  various  reconstruction 
(interpolation)  methods.  Recently,  it  was  suggested  to  min¬ 
imize  the  Sobolev  norm  of  the  reconstruction  error  instead 
of  its  L2  norm;  this  approach  was  then  applied  to  image 
reconstruction  from  singular  points  in  a  Gaussian  scale 
space  [6]. 

However,  within  the  context  of  image  processing  applica¬ 
tions,  the  interpolation  stage  yields  a  discrete-domain  rather 
than  a  continuous-domain  signal  and  no  L2  or  Sobolev 
measures,  but  an  I2  measure,  is  to  be  considered  instead. 

In  this  work,  an  £2  interpolation  method  is  derived  and 
introduced  for  minimizing  the  maximum  possible  approxi¬ 


We  consider  one-dimensional  Sobolev  spaces  defined  over  a 
finite  open  support  D  =  (—7 r,  7 r).  A  Sobolev  space  of  order 
p  is  denoted  by  (^)  and  it  consists  of  all  finite  energy 
functions  for  which  their  first  p  derivatives  are  of  finite 
energy  as  well  [9].  We  adopt  the  following  inner  product 


(x>y)^(fi)  = 


t(xWJW 

n= 0 


l2(q)' 


(2) 


A  Sobolev  space  is  a  reproducing  kernel  Hilbert  space, 
suggesting  an  orthogonal  projection  interpretation  for  the 
ideal  sampling  process. 

Lemma  1:  The  reproducing  kernel  of  H2{fl)  (p  =  1)  is 
given  by 


K(s,t) 


cosh  (| s  —  t\  —  7r) 
2  sinh  (77) 


(3) 


where  s,t  E  Q. 

Proof:  Let  x  E  H2{fl)  be  an  arbitrary  function.  Then, 
it  can  be  expressed  by 


x(s)  =  XA  'eJns’  (4) 

n 

where  equality  holds  point-wise.  The  sampled  value  x(t) 
is  a  linear  bounded  functional  and  by  Riesz  representation 
theorem  can  be  expressed  by  means  of  an  inner  product 
calculation 


x(t)  =  {x(s),K(s,t))H2{ny  (5) 
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Fixing  the  parameter  t,  one  can  express  K(s,t)  by  means 
of  its  Fourier  coefficients  b ,  yielding 


an  •  e* 


xco  =  X 

n 

=  f  x(s)  •  K(s,t)  ds  +  f  x' (s)  •  K'(s,t) 

Jn  Jn 

=  y^nbmjl  +  n-m)  f  ej^n~m^s  ds 

n,m  ^ 

—  27T  ^  ^  bfi  ( 1  H~  Tl  ). 


ds 


(6) 


x  is  arbitrary  yielding  bn  =  (27r)  1  •  e  Jnt/(1  +  n2).  That 
is, 

1  pjn(s-t) 

k^s^  =  ^zY,  ^  •  (?) 


1  +  n2 


Utilizing  the  Fourier  transform  relation  >  2/(l+o;2) 

and  the  aliasing  effect  occurring  in  the  time  domain,  one  can 
show  the  equivalence  of  this  last  expression  with  (3).  □ 

Corollary  1:  The  reproducing  kernel  of  H\  (Q)  is  implic¬ 
itly  given  by 


K(s,t)  = 


1  ^  eMs-t) 

27 r  1  +  n2  +  . . .  +  * 

n 


(8) 


Let  A  =  {£n}n=0  Ar_1  be  a  finite  set  of  sampling 
points.  It  then  follows  that  the  corresponding  sampling 
kernels  {K(s,tn)}n=0  N_1  constitute  a  Riesz  basis  for 

their  span  (the  values  of  a  Sobolev  function  at  distinct  points 
are  linearly  independent) 


S  =  Span{K(s,tn)}^-01.  (9) 

This  sampling  space  is  a  subspace  of  H^iCl).  The  Gram 
matrix  of  these  kernels  can  be  shown  to  comply  with 


G(ra,  n)  =  K(tm ,  tn)  m,  n  =  0, . . . ,  TV  —  1.  (10) 


The  orthogonal  projection  of  x  onto  the  sampling  space  is 
given  by 

N-l 

PsX  =y2an-  (11) 

n=0 

where  a  =  G~1c  and  c  denotes  the  ideal  samples  of  x 
according  to  A.  The  unknown  portion  of  x  that  is  not 
captured  by  the  sampling  process  is  P$^x  =  x  —  P^x. 


III.  INTERPOLATION 

Interpolation  is  the  task  of  evaluating  a  continuous-domain 
function  at  predefined  points  while  having  its  sampled 
version  at  other  points  as  the  only  available  data.  A  com¬ 
mon  approach  to  this  task  relies  on  kernels  that  have 
attractive  properties  in  terms  of  approximation  order,  of 
proportional  constant  and  of  minimal  support  (1).  In  this 
work,  however,  a  different  approach  is  taken  in  which  an  £2 
rather  than  an  L2  measure  is  considered.  Specifically,  the 
interpolated  values  alone  are  compared  with  the  true  values 
while  no  continuous-domain  approximation  is  considered. 
The  motivation  for  such  an  approach  stems  from  image 
processing  applications  for  which  the  interpolation  stage 
yields  a  discrete-domain  rather  than  a  continuous-domain 
signal. 


Theorem  1:  Let  A  =  {tn}n=0  N_1  be  a  sampling  grid 
and  let  c  be  the  ideal  samples  of  x  G  H\  (U)  .  Given  r  ^  A, 
the  solution  of 


argmin  max  |x(r)  —  x(r) 

c,||x||Hp(n)<L 


is  given  by 


x(t)  =  ctG  16, 


(12) 

(13) 


where  G  is  given  by  (10),  bn  =  K(tn,r )  and  K  is  given 
by  (8).  L  is  an  arbitrary  constant  that  complies  with  L2  > 

cTG-1c. 

Proof:  Recalling  Lemma  1,  the  value  x(r)  is  given 
by 

x(r)  =  (x(s),  K(s,t)\  .  (14) 

In  addition,  we  recall  that  A  defines  a  sampling  space  S. 
Therefore,  evaluating  x(r)  is  equivalent  to  the  approxima¬ 
tion  of  (14)  while  having  P^x  as  the  only  available  data. 
This  signal  representation  interpretation  for  the  interpolation 
problem  can  be  applied  now  to  Theorem  1  of  [10].  That  is, 
the  minimax  solution  of  (12)  is  given  by 

x(t)  =  /psx(s),  PsK(s,t)\  .  (15) 

The  sampling  functions  {AT(«s,  tn)}  constitute  a  Riesz  basis 
for  S  and  (13)  follows  accordingly.  The  constant  L,  upper- 
bounding  the  Sobolev  norm  of  x,  is  required  for  defining  a 
robust  minimax  objective  function  [11].  □ 

Theorem  1  describes  a  minimax  approach  to  interpola¬ 
tion.  Given  a  point  to  be  evaluated,  it  defines  an  interpolat¬ 
ing  kernel  to  be  applied  to  the  sampled  data 

k  =  G~1b.  (16) 


The  support  of  this  kernel  may  be  as  large  as  the  size  of 
the  sampled  data;  that  is,  every  sample  value  is  significant 
to  the  interpolation  task.  Nevertheless,  such  kernels  have  in 
practice  a  relatively  small  support  (Fig.  1);  the  reason  for 
that  resides  in  the  structure  of  G-1  which  has  its  significant 
values  located  near  the  main  diagonal.  As  can  be  seen  from 
Fig.  1,  the  support  of  k  increases  as  the  Sobolev  order 
becomes  larger.  Also,  the  ensued  minimax  kernel  for  the 
case  of  p  =  1  corresponds  to  the  average  of  only  the  two 
adjacent  samples  (provided  the  interpolation  point  is  located 
at  the  middle).  The  minimax  solution  of  Theorem  1  is  also 
element-wise  optimal  [11],  implying  that  simultaneously 
designing  several  interpolation  kernels  can  be  performed 
separately,  one  kernel  at  a  time.  Theorem  1  also  gives  rise 
to  an  interpolating  function.  Let  tn  be  the  n-th  sampling 
point.  Then, 

N-l 

4>n(t)  =  X  Gn]n  ■  K(tm,t).  (17) 

m= 0 

This  interpolating  function  is  applied  to  the  sample  value 
originating  from  tn  and  it  consists  of  a  linear  combination 
of  the  sampling  kernels,  i.e.  G  S.  Every  sample  point 
has  its  own  interpolating  function;  if  A  is  uniform,  then 
the  ensued  interpolation  functions  resemble  each  other  and 
posses  a  cyclic  shift-invariant  structure.  Fig.  2  compares 
the  proposed  interpolating  function  with  its  B-spline  coun¬ 
terpart  for  several  Sobolev  orders.  Following  Fig.  2(a),  the 
minimax  approach  yields  an  interpolating  function  that  is 
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(a)  Sobolev  order  p  =  1. 


Fig.  1.  Minimax  interpolation  kernels  for  Sobolev  spaces.  Here,  the 
support  of  the  functions  is  Q  =  (— 7r,  7t),  the  sampling  interval  is  T  =  0.1 
and  the  interpolating  point  is  r  =  0.05.  The  interpolation  kernel  k  is  given 
by  (16)  and  it  is  depicted  here  for  Sobolev  orders  of  pm  1,  2,  3. 


very  similar  to  the  hat  B- spline  function.  However,  these 
two  functions  are  different  as  can  be  observed  from  the 
derivative  of  (j)n  depicted  therein  as  well.  Furthermore,  the 
minimax  approach  yields  finite  support  functions  regardless 
of  the  regularity  constraint  (i.e.,  Sobolev  order)  which  is  in 
contrast  to  B- spline  interpolating  functions.  An  upper  bound 
on  the  approximation  error  is  derived  next. 

Theorem  2:  Let  A  =  {tn}  and  T  =  {rm}  be  a  sampling 
and  an  interpolation  grid,  respectively;  and  let  c  G  RN  and 
d  G  Mm  be  the  ideal  samples  of  x  G  over  A  and  T, 

respectively.  Then, 


d  —  d 


2 

£2 


<B- 


(18) 


where  d  is  the  minimax  approximation  for  d  (Theorem  1), 
G  is  given  by  (10)  and  B  is  the  largest  eigenvalue  of  the 
matrix 


H(k,l)  = 

JV-l 

=  K(rk,ri)-  ^  K(tm,rk)  •  G~1(m,n)  •  K(tn,ri) 

m,n=0 

k,l  =  0...M-  1.  (19) 

Proof:  We  consider  first  a  single  interpolation  point 
and  identify  the  worst-case  input  signal  possible.  Clearly, 
this  signal  satisfies 


x  =  P5x  +  P5ix,  (20) 

where  P^x  is  a  known  continuous-domain  signal  and  is 
given  by  (11).  The  approximation  error  is  given  by 


d  —  d 

= 

(Ps±x,Ps±K(-,r))Hpm 

< 

IP5-lx||hp(0)  •  \\Ps±K(-,t 

where  r  is  the  interpolating  point.  Now, 


llP5-Lx||^p(n) 

—  IMInftn) 

1 1  AT(-,  T-)  1 1  jyP  (fi) 

=  K(  0,0) 

II  PsK^r)\\lm 

=  E K^m 

m,n 

T  sy—1 

c  G  c 


1  Him 


(21) 


,t)  •  G  1(m, n)  ■  K(tn,r), 


(b)  Sobolev  order  p  =  2. 

Fig.  2.  Minimax  interpolating  functions  for  Sobolev  spaces.  Here,  the 
support  of  the  functions  is  Q  =  (— 7r,  7t)  and  the  sampling  interval  is 
T  =  0.5.  The  interpolating  function  is  given  by  (17).  Shown  here  is  a 
comparison  between  the  proposed  interpolating  function  and  the  B -spline 
interpolating  function.  Here,  p  denotes  the  Sobolev  order  of  the  input 
signal. 


and  the  upper  bound  for  this  single  interpolation  point 
follows  immediately  for  the  choice  of 

b  =  ||p5.ir(.,T)||2H|(n) 

=  t)IIh| (fi)  —  \\PsK(’,  t)\\Hp(q)  (22) 

This  upper  bound  is  tight  and  is  achieved  by  signals  of  the 
form 

x  =  P5X  +  a  •  Ps±Kf,r),  (23) 


where  a  is  a  scalar  that  determines  the  Sobolev  norm 
of  x.  When  considering  several  interpolating  points,  the 
aproximation  error  is  given  by 


d-d 


2 

£ 2 


< 


M—l 

E  \(Ps^,PS±K(;Tm))\2 

m= 0 


B-PVxfl 


2 

H%(Q)  » 


(24) 


where  B  is  the  upper  frame  bound  of  the  functions 
{PS±K(-,  Tm)}m=0  M_1  and  is  determined  by  the  largest 
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Fig.  3.  An  upper  bound  on  the  interpolation  error,  given  by  Theorem 
2.  Shown  here  is  B  for  interpolating  by  a  factor  of  two  from  a  uniform 
sampling  grid.  As  the  sampling  interval  shortens,  this  portion  of  the  upper 
bound  tends  to  lower  values.  Here,  p  denotes  the  Sobolev  order  of  the 
input  signal. 


T,  sampling  interval 

Fig.  4.  An  upper  bound  on  the  interpolation  error,  given  by  Theorem 
2.  Shown  here  is  ||Psx||2  =  cTG~1c  for  the  case  of  x  =  (27t)-1/2 
and  for  a  uniform  sampling  grid.  As  the  sampling  interval  shortens,  this 
portion  of  the  upper  bound  approaches  the  value  of  \\x\\ Hv^  =  1,  where 
p  denotes  the  Sobolev  order  of  the  input  signal. 


eigenvalue  of  their  Gram  matrix 

H(k,l)  =  (Ps^Ki-^Ps^K^n))^ 

=  (K(-,Tk),Ps±K(-,Tl))Hp(n) 

=  (K(;  Tk),  K(;  n)  -  PSK( n))H^25) 

It  then  follows  that  the  worst-case  input  signal  achieving 
this  upper  bound  is  of  the  form 

M—l 

X  =  P5X  +  ^2  Pm-  Ps±K(-,  Tm),  (26) 

m=0 


where  /3  is  the  eigenvector  corresponding  to  B  that  also 
guarantees  the  Sobolev  norm  of  x.  □ 

As  the  sampling  interval  shortens,  more  information  on  x 
is  available  and  one  may  expect  the  upper  bound  of  Theorem 
2  to  become  smaller.  This  characteristics  is  manifested  in 
both  B  and  cTG_1c.  Fig.  3  depicts  B  as  a  function  of  the 
sampling  interval  for  a  uniform  grid  and  for  an  interpolation 
by  a  factor  of  two.  Fig.  4  depicts  ||Psx||2  =  cTG~1c  in 
a  similar  manner  for  the  case  of  x  =  (27 r)-1/2;  as  the 
sampling  interval  shortens,  this  portion  of  the  upper  bound 
approaches  the  value  of  ||#||#P(o)  =  1. 

The  generalization  to  images  is  carried  out  by  considering 
2D  Sobolev  functions.  Let  v  =  (^1,^2)  be  a  tuple  of 
nonnegative  integers  where  \v\  =  v%  +  and  let 


Dv 


dUl  dv 2 
da  df3 


(27) 


We  consider  the  two-dimensional  Sobolev  space  of  order 
p  >  2  defined  over  a  finite  open  support  D  =  (— 7r,7r)  x 
(—7 r,  7 r)  (unlike  the  ID  case,  a  2D  Sobolev  space  of  order 
p—  1  is  not  a  reproducing  kernel  Hilbert  space).  This  space 
consists  of  functions  x(a,/?)  that  satisfy  D^x  E  £2(0)  for 
all  possible  v  satisfying  \v\  <  p.  The  corresponding  inner 
product  is  given  by  [9] 


(x-y  >Hf(n)  =  53  <zrx,zry)i2(fi).  (28) 

{u:  0<\u\<p} 


This  inner  product  gives  rise  to  the  following  reproducing 
kernel 


K(a,p,x,y )  =  ^2  53 


ejn(a-x)+jm(p—y) 
^{^:0<| u\<p} 


(29) 


which  is  not  separable  in  a  and  f3.  For  example,  a  Sobolev 
order  of  p  =  2  yields 

l  ^  ejn(cx-x)+jm((3-y) 

(  ’  ^  47 r2  1  +  n2  +  m2  +  n2m2  +  n4  +  m4 

n,m 

(30) 

Next,  the  minimax  approach  is  compared  with  the  cubic 
B-spline  approach.  A  uniform  sampling  grid  is  assumed 
and  an  interpolation  by  a  factor  of  three  is  examined.  Fig. 
5  depicts  a  cervical  Pap  smear  image.  This  image  was 
downsampled  by  a  factor  of  three  and  was  then  interpolated 
accordingly.  The  minimax  interpolation  approach,  shown 
in  Fig.  6,  yields  SNR=30[dB]  while  the  well  known  cu¬ 
bic  B-spline  method  yields  only  29.7[dB].  Table  I  further 
compares  the  minimax  and  the  cubic  B-spline  interpolation 
error  for  several  additional  images.  The  minimax  method 
was  designed  for  minimizing  the  £2  error  and  it  is  shown 
to  yield  better  results  in  terms  of  SNR. 


IV.  CONCLUSIONS 

An  ^2  approach  to  image  reconstruction  has  been  intro¬ 
duced.  The  continuous-domain  images  are  assumed  to  be¬ 
long  to  a  reproducing  kernel  Hilbert  space  and  the  sampling 
process  is  shown  to  correspond  to  an  appropriate  orthogonal 
projection.  It  has  been  also  shown  that  the  values  at  the 
interpolating  grid  correspond  to  a  set  of  inner  product 
calculations,  and  that  this  signal  representation  interpre¬ 
tation  can  be  utilized  to  derive  a  minimax  solution  for 
the  corresponding  1 2  approximation  problem.  To  provide  a 
quantitative  measure  for  the  efficiency  of  the  new  method, 
a  tight  upper  bound  on  the  approximation  error  has  been 
derived  and  demonstrated.  Numerical  examples  show  that 
the  proposed  interpolation  method  yields  better  results  in 
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TABLE  I 

A  COMPARISON  OF  INTERPOLATION  METHODS  BY  A  FACTOR  OF  THREE 


Fig.  5.  A  cervical  Pap  smear  image. 


Fig.  6.  Interpolation  by  a  factor  of  three.  The  proposed  approach,  shown 
here  for  p  =  3,  yields  SNR=30[dB]  while  the  cubic  B-spline  interpolation 
method  yields  only  29.7[dB]. 


terms  of  SNR  compared  to  presently  available  techniques, 
and  could  be  useful  in  various  practical  applications. 
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